#include <stdafx.h>
#include ".\floatvolgrid.h"
#include "stdio.h"

FloatVolGridRaw::FloatVolGridRaw(void)
{
	m_bInterpVoxelsLiterally = false;
	m_nSize[0] = m_nSize[1] = m_nSize[2] = 0;
	m_pData = NULL;
	m_bValid = false;
	m_nNumValidSamples = 0;
	m_lfExpSqrd = 0.0;
	m_lfMax = m_lfMin = m_lfMean = m_lfStdDev = 0.0;
	m_nWorkingChannel = 0;
	m_nNumChannels = 0;
	m_bChannelFirst = false;
}

FloatVolGridRaw::~FloatVolGridRaw(void)
{
	if(m_pData)
	{
		delete [] m_pData;
	}
	m_pData = NULL;
}

int FloatVolGridRaw::GetDataBounds(int nAxis) const
{
	return m_nSize[nAxis];
}

mdsVector<2, double> FloatVolGridRaw::GetSpaceBounds(int nAxis) const
{
	mdsVector<2, double> vecResult;
	vecResult[1] = (double)m_nSize[nAxis];
	return vecResult;
}

float FloatVolGridRaw::GetValAt(int nX, int nY, int nZ) const
{
	return m_pData[GetIdx(nX, nY, nZ, m_nWorkingChannel)];
}

float FloatVolGridRaw::GetValAt(double lfX, double lfY, double lfZ) const
{
	// temp code
	return GetValAt((int)lfX, (int)lfY, (int)lfZ);
}

void FloatVolGridRaw::SetValAt(int nX, int nY, int nZ, float fNewVal)
{
	m_pData[GetIdx(nX, nY, nZ, m_nWorkingChannel)] = fNewVal;
}

mdsVector<3, double> FloatVolGridRaw::EvalGradientAtDataPoint(int nX, int nY, int nZ) const
{
	mdsVector<3, double> vecResult(false);

	if(!m_bValid)
	{
		return vecResult;
	}
	// branch prediction should handle these "if"s appropriately
	if(nX == m_nSize[0]-1) { --nX; }
	if(nY == m_nSize[1]-1) { --nY; }
	if(nZ == m_nSize[2]-1) { --nZ; }
	if(nX < 1) { nX = 1; }
	if(nY < 1) { nY = 1; }
	if(nZ < 1) { nZ = 1; }
	// end bounds checking

	int nIdxBase = GetIdx(nX, nY, nZ, m_nWorkingChannel);
	int nStep = 1;
	if(!m_bChannelFirst)
	{
		nStep = m_nNumChannels;
	}
	vecResult[0] = (double)m_pData[nIdxBase+nStep]-(double)m_pData[nIdxBase-nStep];
	nStep *= m_nSize[0];
	vecResult[1] = (double)m_pData[nIdxBase+nStep]-(double)m_pData[nIdxBase-nStep];
	nStep *= m_nSize[1];
	vecResult[2] = (double)m_pData[nIdxBase+nStep]-(double)m_pData[nIdxBase-nStep];
	// vecResult *= 0.5*GenericMath::oneOverSqrt3;
	return vecResult;
}

mdsVector<3, double> FloatVolGridRaw::EvalGradientAt(double lfX, double lfY, double lfZ) const
{
	// tmp code
	return EvalGradientAtDataPoint((int)lfX, (int)lfY, (int)lfZ);
}

static double DoIntegral(SplineEditCurve * arrSplines, double * arrKnots, int nPts,
						 float * arrPercentiles, int nNumPercentiles,
		double lfMinVal, double lfInvScale)
{
	double lfResult = 0.0;

	static double lfOneThird = (1.0/3.0);

	double lfProbScale = 1.0/nNumPercentiles;

	int iNextPtSpline = 0;
	int iNextPtPercent = 1;

	

	

	double lfXDist;

	double lfNextPt = 0.0;
	double lfCurrPt = 0.0;

	double lfTmp = lfInvScale*(arrPercentiles[0]-lfMinVal);
	while(arrKnots[iNextPtSpline] < lfTmp && iNextPtSpline < nPts)
	{
		lfCurrPt = arrKnots[iNextPtSpline];
		++iNextPtSpline;
	}
	if(lfTmp > lfCurrPt)
	{
		lfCurrPt = lfTmp;
	}


	

	while(iNextPtSpline < 1 && iNextPtPercent < nNumPercentiles)
	{
		while(lfInvScale*(arrPercentiles[iNextPtPercent]-lfMinVal) < 0.0 
			&& iNextPtPercent < nNumPercentiles)
		{
			++iNextPtPercent;
		}
		if(iNextPtPercent >= nNumPercentiles)
		{
			return 0.0;
		}

		lfNextPt = lfInvScale*(arrPercentiles[iNextPtPercent]-lfMinVal);

		if(arrKnots[iNextPtSpline] < lfNextPt)
		{
			lfNextPt = arrKnots[iNextPtSpline];
		}

		

		// now do stuff on appropriate interval

		lfXDist = lfNextPt-lfCurrPt;

		double lfScaleIvalBy = arrPercentiles[iNextPtPercent] - 
			arrPercentiles[iNextPtPercent-1];
		lfScaleIvalBy *= lfInvScale; 

		lfResult += (arrSplines[iNextPtSpline].IntegrateIval(lfCurrPt, lfNextPt))*lfProbScale/lfScaleIvalBy;
			

		lfCurrPt = lfNextPt;

		if(arrKnots[iNextPtSpline] < lfInvScale*(arrPercentiles[iNextPtPercent]-lfMinVal))
		{
			++iNextPtSpline;
		}
		else
		{
			++iNextPtPercent;
		}
	}	

	while(iNextPtSpline < nPts && iNextPtPercent < nNumPercentiles)
	{
		double lfNextPt = lfInvScale*(arrPercentiles[iNextPtPercent]-lfMinVal);

		if(arrKnots[iNextPtSpline] < lfNextPt)
		{
			lfNextPt = arrKnots[iNextPtSpline];
		}

		// now do stuff on appropriate interval

		lfXDist = lfNextPt-lfCurrPt;

		double lfScaleIvalBy = arrPercentiles[iNextPtPercent] - 
			arrPercentiles[iNextPtPercent-1];
		lfScaleIvalBy *= lfInvScale; 

		lfResult += (arrSplines[iNextPtSpline].IntegrateIval(lfCurrPt, lfNextPt))*lfProbScale/lfScaleIvalBy;

		lfCurrPt = lfNextPt;

		if(arrKnots[iNextPtSpline] < lfInvScale*(arrPercentiles[iNextPtPercent]-lfMinVal))
		{
			++iNextPtSpline;
		}
		else
		{
			++iNextPtPercent;
		}
	}	
	//if(iNextPtSpline



	return lfResult;
}

bool FloatVolGridRaw::DoProbXferFunction(FloatVolGridRaw & src, 
		SplineEditCurve redPolys[], double redKnots[], int nRedPts,
		SplineEditCurve greenPolys[], double greenKnots[], int nGreenPts,
		SplineEditCurve bluePolys[], double blueKnots[], int nBluePts, 
		SplineEditCurve alphaPolys[], double alphaKnots[], int nAlphaPts,
		double lfMinVal, double lfScale)
{
	if(src.m_nSize[0] != m_nSize[0] || 
		src.m_nSize[1] != m_nSize[1] || 
		src.m_nSize[2] != m_nSize[2])
	{
		return false;
	}

	double lfInvScale = 1.0/lfScale;

	if(m_bChannelFirst || src.m_bChannelFirst)
	{
		return false;
	}

	double lfProbScale = 1.0/src.m_nNumChannels;

	for(int k = 0; k < src.m_nSize[2]; ++k)
	{
		for(int j = 0; j < src.m_nSize[1]; ++j)
		{
			for(int i = 0; i < src.m_nSize[0]; ++i)
			{
				float * arrFPercenties = &src.m_pData[src.GetIdx(i,j,k, 0)];


				float * arrFColors = &m_pData[GetIdx(i,j,k,0)];

				if(arrFPercenties[0] < -100)
				{
					int nCnt = 0;
					int iMin = (i-1)*(i >= 1);
					int iMax = (i+2-m_nSize[0])*((i+2) <= m_nSize[0]) + m_nSize[0];

					int jMin = (j-1)*(j >= 1);
					int jMax = (j+2-m_nSize[1])*((j+2) <= m_nSize[1]) + m_nSize[1];

					int kMin = (k-1)*(k >= 1);
					int kMax = (k+2-m_nSize[2])*((k+2) <= m_nSize[2]) + m_nSize[2];

					int nCntPossible = (2+(i < 1)+(i+2>m_nSize[0]))*
						(2+(j < 1)+(j+2>m_nSize[1]))*
						(2+(k < 1)+(k+2>m_nSize[2]));

					for(int kTmp = kMin; kTmp < kMax; ++kTmp)
					{
						for(int jTmp = jMin; jTmp < jMax; ++jTmp)
						{
							for(int iTmp = iMin; iTmp < iMax; ++iTmp)
							{
								nCnt += (src.m_pData[GetIdx(iTmp, jTmp, kTmp, 0)] < -100);
							}
						}
					}
					arrFColors[3] = 0.8+0.2*nCnt/((double)nCntPossible);
					arrFColors[2] = 0.75;
					arrFColors[1] = 0.75;
					arrFColors[0] = 0.75;
				}
				else
				{

					arrFColors[0] = DoIntegral(redPolys, redKnots, nRedPts, arrFPercenties, 
						src.m_nNumChannels, lfMinVal, lfInvScale);

					arrFColors[1] = DoIntegral(greenPolys, greenKnots, nGreenPts, arrFPercenties, 
						src.m_nNumChannels, lfMinVal, lfInvScale);

					arrFColors[2] = DoIntegral(bluePolys, blueKnots, nBluePts, arrFPercenties, 
						src.m_nNumChannels, lfMinVal, lfInvScale);

					arrFColors[3] = 2.0*DoIntegral(alphaPolys, alphaKnots, nAlphaPts, arrFPercenties, 
						src.m_nNumChannels, lfMinVal, lfInvScale);
				}
			}
		}
	}

	m_bValid = true;
	return true;
}

bool FloatVolGridRaw::SetSize(int nX, int nY, int nZ, int nChannels)
{
	m_nSize[0] = nX;
	m_nSize[1] = nY;
	m_nSize[2] = nZ;
	int nBytesPerDatEl = sizeof(float)*nChannels;
	int nFloatsPerVoxel = nChannels;
	m_nNumChannels = nFloatsPerVoxel;
	m_pData = new float[m_nSize[0]*m_nSize[1]*m_nSize[2]*m_nNumChannels];

	if(!m_pData)
	{
		m_bValid = false;
		return false;
	}

	return true;
}

bool FloatVolGridRaw::ConvertFromRawVolGrid(const RawVolGrid & rawVolGrid, bool bIsBigEndian, bool bFlip)
{
	m_nSize[0] = rawVolGrid.GetSizeX();
	m_nSize[1] = rawVolGrid.GetSizeY();
	m_nSize[2] = rawVolGrid.GetSizeZ();
	m_nNumChannels = rawVolGrid.GetVoxelSize();
	m_nNumChannels /= 4;
	m_pData = new float[m_nSize[0]*m_nSize[1]*m_nSize[2]*m_nNumChannels];

	if(!m_pData)
	{
		m_bValid = false;
		return false;
	}

	double lfSumValidVals = 0.0;
	double lfSumSqrdValidVals = 0.0;
	int nCntValidSamples = 0;
	for(int k = 0; k < m_nSize[2]; ++k)
	{
		for(int j = 0; j < m_nSize[1]; ++j)
		{
			int iStart = 0;
			int iStop = m_nSize[0]-1;
			int iStep = 1;
			if(bFlip)
			{
				iStep = -1;
				iStart = iStop;
				iStop = 0;
			}
			int iSrc = 0;
			for(int i = iStart; i != iStop; i += iStep) // dangerous use of "!=" rather then <
			{
				for(int nChannel = 0; nChannel < m_nNumChannels; ++nChannel)
				{
					unsigned char pcEndianReverse[4];
					const unsigned char * pSrc = rawVolGrid.GetAt(iSrc,j,k)+4*nChannel;
					if(bIsBigEndian)
					{
						pcEndianReverse[3] = pSrc[0];
						pcEndianReverse[2] = pSrc[1];
						pcEndianReverse[1] = pSrc[2];
						pcEndianReverse[0] = pSrc[3];
					}
					else
					{
						pcEndianReverse[0] = pSrc[0];
						pcEndianReverse[1] = pSrc[1];
						pcEndianReverse[2] = pSrc[2];
						pcEndianReverse[3] = pSrc[3];
					}
					float * pfVal = (float *)((unsigned char *)pcEndianReverse);

					float fTest = 5.2e2f;
					if(*pfVal > -100.0) // *pfVal < 1.0e35f && *pfVal > -1.0e35f)
					{
						if(i == 0 && j == 0 && k == 0)
						{
							m_lfMax = m_lfMin = *pfVal;
						}

						if(*pfVal > m_lfMax)
						{
							m_lfMax = *pfVal;
						}
						if(*pfVal < m_lfMin)
						{
							m_lfMin = *pfVal;
						}

						lfSumValidVals += (double)(*pfVal);
						lfSumSqrdValidVals += ((double)(*pfVal))*((double)(*pfVal));
						++nCntValidSamples;
					}
					
					SetWorkingChannel(nChannel);
					SetValAt(i,j,k, *pfVal);
					
				}
				++iSrc;
			}
		}
	}
	SetWorkingChannel(0);
	m_bValid = (m_nSize[0] > 1) && (m_nSize[1] > 1) && (m_nSize[2] > 1) && (nCntValidSamples > 0);
	m_lfMean = lfSumValidVals/nCntValidSamples;
	m_lfExpSqrd = lfSumSqrdValidVals/nCntValidSamples;
	m_lfStdDev = m_lfExpSqrd - m_lfMean*m_lfMean;
	m_lfStdDev = sqrt(fabs(m_lfStdDev));
	m_nNumValidSamples = nCntValidSamples;
	return true;
}


///////////////////////////////////////////////////////////////////////////////////////////
// FloatVolGridDownSampled
FloatVolGridDownSampled::FloatVolGridDownSampled(void)
{
	m_nSize[0] = m_nSize[1] = m_nSize[2] = 0;
	m_nWorkingChannel = 0;
	m_pData = NULL;
	m_bValid = false;
	m_nNumValidSamples = 0;
	m_lfExpSqrd = 0.0;
	m_lfMax = m_lfMin = m_lfMean = m_lfStdDev = 0.0;
	
}

FloatVolGridDownSampled::~FloatVolGridDownSampled(void)
{
	if(m_pData)
	{
		delete [] m_pData;
	}
	m_pData = NULL;
}

int FloatVolGridDownSampled::GetDataBounds(int nAxis) const
{
	return 2*m_nSize[nAxis];
}

mdsVector<2, double> FloatVolGridDownSampled::GetSpaceBounds(int nAxis) const
{
	mdsVector<2, double> vecResult;
	vecResult[1] = 2.0*(double)m_nSize[nAxis];
	return vecResult;
}

float FloatVolGridDownSampled::GetValAt(int nX, int nY, int nZ) const
{
	return m_pData[GetIdx(nX>>1, nY>>1, nZ>>1)];
	//int nXOneUp = (nX+1)>>1;
	//int nYOneUp = (nY+1)>>1;
	//int nZOneUp = (nZ+1)>>1;
	//nX >>= 1;
	//nY >>= 1;
	//nZ >>= 1;
	//nXOneUp -= (nXOneUp >= m_nSize[0]);
	//nYOneUp -= (nYOneUp >= m_nSize[1]);
	//nZOneUp -= (nZOneUp >= m_nSize[2]);
	//// int nIdxLookup = GetIdx(
	//return 0.125f*(
	//	m_pData[GetIdx(nX     , nY     , nZ     )] +
	//	m_pData[GetIdx(nX     , nY     , nZOneUp)] +
	//	m_pData[GetIdx(nX     , nYOneUp, nZ     )] +
	//	m_pData[GetIdx(nX     , nYOneUp, nZOneUp)] +
	//	m_pData[GetIdx(nXOneUp, nY     , nZ     )] +
	//	m_pData[GetIdx(nXOneUp, nY     , nZOneUp)] +
	//	m_pData[GetIdx(nXOneUp, nYOneUp, nZ     )] +
	//	m_pData[GetIdx(nXOneUp, nYOneUp, nZOneUp)]
	//	);
}

float FloatVolGridDownSampled::GetValAt(double lfX, double lfY, double lfZ) const
{
	// temp code
	return GetValAt((int)lfX, (int)lfY, (int)lfZ);
}

void FloatVolGridDownSampled::SetValAt(int nX, int nY, int nZ, float fNewVal)
{
	m_pData[GetIdx(nX>>1, nY>>1, nZ>>1)] = fNewVal;
}

mdsVector<3, double> FloatVolGridDownSampled::EvalGradientAtDataPoint(int nX, int nY, int nZ) const
{
	mdsVector<3, double> vecResult(false);

	if(!m_bValid)
	{
		return vecResult;
	}
	nX >>=1;
	nY >>=1;
	nZ >>=1;
	// branch prediction should handle these "if"s appropriately
	if(nX == m_nSize[0]-1) { --nX; }
	if(nY == m_nSize[1]-1) { --nY; }
	if(nZ == m_nSize[2]-1) { --nZ; }
	if(nX < 1) { nX = 1; }
	if(nY < 1) { nY = 1; }
	if(nZ < 1) { nZ = 1; }
	// end bounds checking

	int nIdxBase = GetIdx(nX, nY, nZ);
	vecResult[0] = (double)m_pData[nIdxBase+1]-(double)m_pData[nIdxBase-1];
	int nStep = m_nSize[0];
	vecResult[1] = (double)m_pData[nIdxBase+nStep]-(double)m_pData[nIdxBase-nStep];
	nStep *= m_nSize[1];
	vecResult[2] = (double)m_pData[nIdxBase+nStep]-(double)m_pData[nIdxBase-nStep];
	vecResult *= 0.5; //*GenericMath::oneOverSqrt3;
	return vecResult;
}

mdsVector<3, double> FloatVolGridDownSampled::EvalGradientAt(double lfX, double lfY, double lfZ) const
{
	// tmp code
	return EvalGradientAtDataPoint((int)lfX, (int)lfY, (int)lfZ);
}

bool FloatVolGridDownSampled::LoadFromPreprocessedFile(const char * fName, int nX, int nY, int nZ)
{
	// return false;
	nX >>= 1; nY >>= 1; nZ >>= 1;
	FILE * fIn = fopen(fName, "rb");
	if(!fIn)
	{
		return false;
	}

	if(m_pData)
	{
		delete [] m_pData;
	}
	m_pData = new float[nX*nY*nZ];

	size_t nBytesRead = fread(m_pData, sizeof(float), nX*nY*nZ, fIn);
	if(nBytesRead < (size_t)(nX*nY*nZ))
	{
		delete [] m_pData;
		m_pData = NULL;
		fclose(fIn);
		return false;
	}

	m_nSize[0] = nX;
	m_nSize[1] = nY;
	m_nSize[2] = nZ;

	if(fread(&m_lfExpSqrd, sizeof(double), 1, fIn) < 1 || 
		fread(&m_lfMean, sizeof(double), 1, fIn) < 1 || 
		fread(&m_nNumValidSamples, sizeof(int), 1, fIn) < 1 || 
		fread(&m_lfMin, sizeof(double), 1, fIn) < 1 || 
		fread(&m_lfMax, sizeof(double), 1, fIn) < 1)
	{
		delete [] m_pData;
		m_pData = NULL;
		fclose(fIn);
		return false;
	}

	//m_lfExpSqrd = m_pData[m_nSize[0]*m_nSize[1]*m_nSize[2]];
	//m_lfMean = m_pData[m_nSize[0]*m_nSize[1]*m_nSize[2]+1];
	//m_nNumValidSamples = (int)m_pData[m_nSize[0]*m_nSize[1]*m_nSize[2]+2];
	//m_lfMin = m_pData[m_nSize[0]*m_nSize[1]*m_nSize[2]+3];
	//m_lfMax = m_pData[m_nSize[0]*m_nSize[1]*m_nSize[2]+4];

	m_lfStdDev = m_lfExpSqrd - m_lfMean*m_lfMean;
	m_lfStdDev = sqrt(fabs(m_lfStdDev));

	fclose(fIn);

	m_bValid = true;

	return true;
}

bool FloatVolGridDownSampled::DumpToPreprocessedFile(const char * fName) const
{
	// return false;
	FILE * fOut = fopen(fName, "wb");
	if(!fOut)
	{
		return false;
	}
	

	size_t numFloatsToWrite = m_nSize[0]*m_nSize[1]*m_nSize[2];

	size_t numFloatsWritten = fwrite(m_pData, sizeof(float), numFloatsToWrite, fOut);

	if(numFloatsWritten < numFloatsToWrite)
	{
		fclose(fOut);
		return false;
	}

	if(fwrite(&m_lfExpSqrd, sizeof(double), 1, fOut) < 1 || 
		fwrite(&m_lfMean, sizeof(double), 1, fOut) < 1 || 
		fwrite(&m_nNumValidSamples, sizeof(int), 1, fOut) < 1 || 
		fwrite(&m_lfMin, sizeof(double), 1, fOut) < 1 || 
		fwrite(&m_lfMax, sizeof(double), 1, fOut) < 1)
	{
		fclose(fOut);
		return false;
	}

	fclose(fOut);
	return true;
}




bool FloatVolGridDownSampled::ConvertFromRawVolGrid(const RawVolGrid & rawVolGrid, bool bIsBigEndian, bool bFlip)
{
	m_nSize[0] = rawVolGrid.GetSizeX();
	m_nSize[1] = rawVolGrid.GetSizeY();
	m_nSize[2] = rawVolGrid.GetSizeZ();
	int nBytesPerDatEl = rawVolGrid.GetVoxelSize();
	int nFloatsPerVoxel = (nBytesPerDatEl/4)+((nBytesPerDatEl%4) != 0);
	m_nSize[0] >>= 1;
	m_nSize[1] >>= 1;
	m_nSize[2] >>= 1;
	m_pData = new float[m_nSize[0]*m_nSize[1]*m_nSize[2]*nFloatsPerVoxel]; // added space for 5 floats for statistics

	if(!m_pData)
	{
		m_bValid = false;
		return false;
	}

	double lfSumValidVals = 0.0;
	double lfSumSqrdValidVals = 0.0;
	int nCntValidSamples = 0;
	for(int k = 0; k < 2*m_nSize[2]; ++k)
	{
		for(int j = 0; j < 2*m_nSize[1]; ++j)
		{
			for(int i = 0; i < 2*m_nSize[0]; ++i)
			{
				unsigned char pcEndianReverse[4];
				const unsigned char * pSrc = rawVolGrid.GetAt(i,j,k);
				if(bIsBigEndian)
				{
					pcEndianReverse[3] = pSrc[0];
					pcEndianReverse[2] = pSrc[1];
					pcEndianReverse[1] = pSrc[2];
					pcEndianReverse[0] = pSrc[3];
				}
				else
				{
					pcEndianReverse[0] = pSrc[0];
					pcEndianReverse[1] = pSrc[1];
					pcEndianReverse[2] = pSrc[2];
					pcEndianReverse[3] = pSrc[3];
				}
				float * pfVal = (float *)((unsigned char *)pcEndianReverse);

				float fTest = 5.2e2f;
				if(*pfVal < 1.0e35f && *pfVal > -1.0e35f)
				{
					if(i == 0 && j == 0 && k == 0)
					{
						m_lfMax = m_lfMin = *pfVal;
					}

					if(*pfVal > m_lfMax)
					{
						m_lfMax = *pfVal;
					}
					if(*pfVal < m_lfMin)
					{
						m_lfMin = *pfVal;
					}

					lfSumValidVals += (double)(*pfVal);
					lfSumSqrdValidVals += ((double)(*pfVal))*((double)(*pfVal));
					++nCntValidSamples;
				}
				
				if( ! ((i%2) || (j%2) || (k%2)))
				{
					SetValAt(i,j,k, 0.125f*(*pfVal));
				}
				else
				{
					m_pData[GetIdx(i>>1, j>>1, k>>1)] += 0.125f*(*pfVal);
				}
			}
		}
	}
	m_bValid = (m_nSize[0] > 1) && (m_nSize[1] > 1) && (m_nSize[2] > 1) && (nCntValidSamples > 0);
	m_lfMean = lfSumValidVals/nCntValidSamples;
	m_lfExpSqrd = lfSumSqrdValidVals/nCntValidSamples;
	m_lfStdDev = m_lfExpSqrd - m_lfMean*m_lfMean;
	m_lfStdDev = sqrt(fabs(m_lfStdDev));
	m_nNumValidSamples = nCntValidSamples;

	// store them here as well, for quick access when reading/writing
	//m_pData[m_nSize[0]*m_nSize[1]*m_nSize[2]] = (float)m_lfExpSqrd;
	//m_pData[m_nSize[0]*m_nSize[1]*m_nSize[2]+1] = (float)m_lfMean;
	//m_pData[m_nSize[0]*m_nSize[1]*m_nSize[2]+2] = (float)m_nNumValidSamples;
	//m_pData[m_nSize[0]*m_nSize[1]*m_nSize[2]+3] = (float)m_lfMin;
	//m_pData[m_nSize[0]*m_nSize[1]*m_nSize[2]+4] = (float)m_lfMax;

	return true;
}

// VolumeViz implementation
#include <gl/gl.h>

viewInfo VolumeViz::m_gDummyViewInfo;

VolumeViz::VolumeViz():m_nRenderMode(1)
{
	m_nStepSize[0] = m_nStepSize[1] = m_nStepSize[2] = 4;
	m_pVolGrid = NULL;
	m_pViewInfo = &VolumeViz::m_gDummyViewInfo;
	m_lfTransferOffset = 0.0;
	m_lfTransferWidthScale = 1.0;
}

VolumeViz::VolumeViz(FloatVolGrid & volGrid):m_nRenderMode(1)
{
	m_nStepSize[0] = m_nStepSize[1] = m_nStepSize[2] = 4;
	SetVolGrid(volGrid);
	m_pViewInfo = &VolumeViz::m_gDummyViewInfo;
	m_lfTransferOffset = 0.0;
	m_lfTransferWidthScale = 1.0;
}

VolumeViz::~VolumeViz()
{
}

void VolumeViz::Draw()
{
	if(!m_pVolGrid)
	{
		return;
	}

	static GLfloat lclLightSpecular[] = { 1.0f, 1.0f, 1.0f, 1.0f };
	static GLfloat lclLightAmbDiffuse[] = { 0.2f, 0.0f, 0.7f, 1.0f };
	glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR,lclLightSpecular);
	glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, lclLightAmbDiffuse);
	glMateriali(GL_FRONT_AND_BACK, GL_SHININESS,128);

	bool bLiteralVoxels = m_pVolGrid->GetVoxelInterpretation();

	glPushMatrix();

	
	glTranslated(-m_pVolGrid->GetDataBounds(0)*0.5, -m_pVolGrid->GetDataBounds(1)*0.5,
		-m_pVolGrid->GetDataBounds(2)*0.5);
	
	mdsVector<3, double> vecXlate;
	vecXlate[0] = -m_pVolGrid->GetDataBounds(0)*0.5;
	vecXlate[1] = -m_pVolGrid->GetDataBounds(1)*0.5;
	vecXlate[2] = -m_pVolGrid->GetDataBounds(2)*0.5;

	glEnable(GL_BLEND);
	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
	
	if(!m_nRenderMode)
	{
//#if !USE_HEX_SPRITES
		glBegin(GL_QUADS);// old quad-based sprites
//#endif
	}

	
	double lfOneOverTwiceStdDev = 0.5/
		(m_lfStdDev);
	// 0.5/(vizGlobals::cloudGrids[vizGlobals::nAnimFrame].GetStdDev());
	
	int xStep = m_nStepSize[0];
	int yStep = m_nStepSize[1];
	int hStep = m_nStepSize[2];
	int nMin[3];
	int nMax[3];
	int nAxis = 0;
	nMin[0] = nMin[1] = nMin[2] = 0;
	for(nAxis = 0; nAxis < 3; ++nAxis)
	{
		nMax[nAxis] = m_pVolGrid->GetDataBounds(nAxis);
	}
	
	int nViewDirQuad = OctUtils::Vec3To3dQuadrant(m_pViewInfo->GetLookDir());
	mdsVector<3, double> vecCurrViewDir = m_pViewInfo->GetLookDir();
	mdsVector<3, double> vecCurrViewPos = m_pViewInfo->GetPosition()-vecXlate;
	mdsVector<3, double> vecUp = m_pViewInfo->GetUpDir();
	mdsVector<3, double> vecLeft = m_pViewInfo->GetLeftDir();

	
	

	//double lfNewDist = fabs(vecCurrViewDir[0]*(nMin[0]-nMax[0])) + 
	//	fabs(vecCurrViewDir[1]*(nMin[1]-nMax[1])) + 
	//	fabs(vecCurrViewDir[2]*(nMin[2]-nMax[2]));
	//double lfNewAlpha = 0.3; //*(exp(100.0/lfNewDist));
	for(int nWhichQuad = 0; nWhichQuad < 8 ;++nWhichQuad)
	{
		int nLclQuad = nViewDirQuad^(OctUtils::QuadTraverseOrder[nWhichQuad]);
		int nStart[3];
		int nStep[3];
		int nStop[3];
		int nCompSignum[3];
		bool bDoThisLoop = true;


	mdsVector<3, double> arrVecHexPoints[7];
	int iHexPoint = 0;
	for(;iHexPoint < 7; ++iHexPoint)
	{
		int iCurrCorner = nLclQuad^(OctUtils::OctCornerGreyCode[iHexPoint+1]);
		// int iCurrCorner = nViewDirQuad^(OctUtils::OctCornerGreyCode[iHexPoint+1]);
		for(nAxis = 0; nAxis < 3; ++nAxis)
		{
			if(iCurrCorner&1)
			{
				arrVecHexPoints[iHexPoint][nAxis] = 0.5*m_nStepSize[nAxis];
			}
			else
			{
				arrVecHexPoints[iHexPoint][nAxis] = -0.5*m_nStepSize[nAxis];
			}
			iCurrCorner >>= 1;
		}
		
		mdsVector<3, double> stepBackBy(m_pViewInfo->GetLookDir());
		double lfStepBackDistance = arrVecHexPoints[iHexPoint].Dot(vecCurrViewDir);
		stepBackBy *= lfStepBackDistance;
		// arrVecHexPoints[iHexPoint] -= stepBackBy;
	}
	// don't worry about face orientation, as we set the normals manually, and as
	// backface culling would be silly in the absence of any actual backfaces

		for(nAxis = 0; nAxis < 3; ++nAxis)
		{
			int mask = 1<<nAxis;
			int nBump = (nLclQuad&mask)>>nAxis;
			nStop[nAxis] = static_cast<int>(vecCurrViewPos[nAxis]);
			if(mask & nLclQuad)
			{
				nStart[nAxis] = nMax[nAxis]-1;
				nStep[nAxis] = -m_nStepSize[nAxis];
				nStop[nAxis] -= (m_nStepSize[nAxis]-nBump);
				
				if(nStop[nAxis] < nMin[nAxis])
				{
					nStop[nAxis] = nMin[nAxis]-1;
				}
				
				if(nStop[nAxis] >= nStart[nAxis])
				{
					bDoThisLoop = false;
				}
				nCompSignum[nAxis] = 1;
				
			}
			else
			{
				nStart[nAxis] = nMin[nAxis];
				nStep[nAxis] = m_nStepSize[nAxis];
				nStop[nAxis] += (m_nStepSize[nAxis]-nBump);
				if(nStop[nAxis] >= nMax[nAxis])
				{
					nStop[nAxis] = nMax[nAxis];
				}
				
				if(nStop[nAxis] <= nStart[nAxis])
				{
					bDoThisLoop = false;
				}
				nCompSignum[nAxis] = -1;
			}
		}
		if(!bDoThisLoop)
		{
			continue;
		}
		
		for(int hVal = nStart[2]; (hVal-nStop[2])*nCompSignum[2] > 0; hVal += nStep[2])
		{
			for(int yVal= nStart[1]; (yVal-nStop[1])*nCompSignum[1] > 0; yVal += nStep[1])
			{
				for(int xVal = nStart[0]; (xVal-nStop[0])*nCompSignum[0] > 0; xVal += nStep[0])
				{
					
					// m_pVolGrid->EvalGradientAt(
					float fCurrVal, fRed, fGreen, fBlue, fAlpha;
					if(bLiteralVoxels)
					{
						// there has got to be a more efficient way....
						m_pVolGrid->SetWorkingChannel(0);
						fRed = 
							m_pVolGrid->GetValAt(xVal, yVal, hVal);
						m_pVolGrid->SetWorkingChannel(1);
						fGreen = 
							m_pVolGrid->GetValAt(xVal, yVal, hVal);
						m_pVolGrid->SetWorkingChannel(2);
						fBlue = 
							m_pVolGrid->GetValAt(xVal, yVal, hVal);

						m_pVolGrid->SetWorkingChannel(3); // for alpha
						fAlpha = 
							m_pVolGrid->GetValAt(xVal, yVal, hVal);
					}
					else
					{
						fCurrVal = 
							m_pVolGrid->GetValAt(xVal, yVal, hVal);

						// MDS_MORE_CODE_HERE
						//if(vizGlobals::renderInfo::bSkipBadValues && fabs(fCurrVal) >= 1.0e33f)
						//{
						//	continue;
						//}

						
						fCurrVal = (float)(
							(fCurrVal+m_lfStdDev -
							m_lfMean)*
							lfOneOverTwiceStdDev);
					}

					

					mdsVector<3, double> gradient =
						m_pVolGrid->EvalGradientAt(xVal, yVal, hVal);


					if(m_nRenderMode)
					{
						glBegin(GL_TRIANGLE_FAN);
					}

					gradient.Normalize();
					gradient *= -1.0;
					::glNormal3d(gradient[0], gradient[1], gradient[2]);

					gradient *= -1.0;
					
					if(bLiteralVoxels)
					{
						glColor4f(fRed, fGreen, fBlue, fAlpha);
					}
					else
					{
						float fClampedOpac = 0.5f+(fCurrVal-0.5f-
							m_lfTransferOffset)*
							m_lfTransferWidthScale;
						if(fClampedOpac > 1.0f)
						{
							fClampedOpac = 1.0f;
						}
						else if(fClampedOpac < 0.0f)
						{
							fClampedOpac = 0.0f;
						}
						else
						{
							fClampedOpac = fClampedOpac*fClampedOpac*(3-2*fClampedOpac); // 3x^2-2x^3
						}
						fRed = fCurrVal; fGreen = fCurrVal; fBlue = 1.0f-fCurrVal;
						fAlpha = 0.3*fClampedOpac;
						glColor4f(fCurrVal, fCurrVal, 1.0f-fCurrVal, 0.3*fClampedOpac); //0.3f*fCurrVal);
					}

					if(m_nRenderMode)
					{

						glVertex3f(
							(GLfloat)(xVal+arrVecHexPoints[6].x()),
							(GLfloat)(yVal+arrVecHexPoints[6].y()),
							(GLfloat)(hVal+arrVecHexPoints[6].z()));
						glColor4f(fRed, fGreen, fBlue, 0.1*fAlpha);
						for(int nHexPoint = 5; nHexPoint >= 0; --nHexPoint)
						{
							glVertex3f(
								(GLfloat)(xVal+arrVecHexPoints[nHexPoint].x()),
								(GLfloat)(yVal+arrVecHexPoints[nHexPoint].y()),
								(GLfloat)(hVal+arrVecHexPoints[nHexPoint].z()));
						}
						glVertex3f(
							(GLfloat)(xVal+arrVecHexPoints[5].x()),
							(GLfloat)(yVal+arrVecHexPoints[5].y()),
							(GLfloat)(hVal+arrVecHexPoints[5].z()));

						glEnd(); // end of new hexagon-based code for sprites
					}
					else
					{
						// // old quad-based sprites
						//
						glVertex3f(
							(GLfloat)(xVal+0.5*xStep*vecUp.x()-0.5*xStep*vecLeft.x()), 
							(GLfloat)(yVal+0.5*yStep*vecUp.y()-0.5*yStep*vecLeft.y()), 
							(GLfloat)(hVal+0.5*hStep*vecUp.z()-0.5*hStep*vecLeft.z())); // up right
						glVertex3f(
							(GLfloat)(xVal+0.5*xStep*vecUp.x()+0.5*xStep*vecLeft.x()), 
							(GLfloat)(yVal+0.5*yStep*vecUp.y()+0.5*yStep*vecLeft.y()), 
							(GLfloat)(hVal+0.5*hStep*vecUp.z()+0.5*hStep*vecLeft.z())); // up left
						glVertex3f(
							(GLfloat)(xVal-0.5*xStep*vecUp.x()+0.5*xStep*vecLeft.x()), 
							(GLfloat)(yVal-0.5*yStep*vecUp.y()+0.5*yStep*vecLeft.y()), 
							(GLfloat)(hVal-0.5*hStep*vecUp.z()+0.5*hStep*vecLeft.z())); // down left
						glVertex3f(
							(GLfloat)(xVal-0.5*xStep*vecUp.x()-0.5*xStep*vecLeft.x()), 
							(GLfloat)(yVal-0.5*yStep*vecUp.y()-0.5*yStep*vecLeft.y()), 
							(GLfloat)(hVal-0.5*hStep*vecUp.z()-0.5*hStep*vecLeft.z())); // down right

						// end of code for old quad-based sprites
					}


				}
			}
		}
	}

	if(!m_nRenderMode)
	{
		glEnd(); // old quad-based sprites
	}

	glPopMatrix();
}